

In [1]:

    
%load_ext load_style
%load_style talk.css

EOF analysis - global hgt500

In statistics and signal processing, the method of empirical orthogonal function (EOF) analysis is a decomposition of a signal or data set in terms of orthogonal basis functions which are determined from the data. It is similar to performing a principal components analysis on the data, except that the EOF method finds both temporal projections and spatial patterns. The term is also interchangeable with the geographically weighted PCAs in geophysics (https://en.wikipedia.org/wiki/Empirical_orthogonal_functions). The spatial patterns are the EOFs, and can be thought of as basis functions in terms of variance. The associated temporal projections are the pricipal components (PCs) and are the temporal coefficients of the EOF patterns.

Python EOF library

The EOF library is used here from Dawson (http://ajdawson.github.io/eofs/). eofs is a Python package for EOF analysis of spatial-temporal data. Using EOFs (empirical orthogonal functions) is a common technique to decompose a signal varying in time and space into a form that is easier to interpret in terms of spatial and temporal variance. Some of the key features of eofs are:

Suitable for large data sets: computationally efficient for the large output data sets of modern climate models.

Transparent handling of missing values: missing values are removed automatically during computations and placed back into output fields.

Automatic metadata: metadata from input fields is used to construct metadata for output fields.

No Compiler required: a fast implementation written in pure Python using the power of numpy, no Fortran or C dependencies.

Data Source

The hgt data is downloaded from https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.pressure.html. We select hgt at 500mb level from 1979 to 2003, and convert them into yearly climatology using CDO.

cdo yearmean -sellevel,500/500, -selyear,1979/2003 hgt.mon.mean.nc hgt500.mon.mean.nc

1. Load basic libraries



In [2]:

    
% matplotlib inline

import numpy as np
from scipy import signal
import numpy.polynomial.polynomial as poly
from netCDF4 import Dataset

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from eofs.standard import Eof

2. Load ght500 data



In [3]:

    
infile = 'data/hgt500.mon.mean.nc'
ncin = Dataset(infile, 'r')
hgt  = ncin.variables['hgt'][:,0,:,:]
lat  = ncin.variables['lat'][:]
lon  = ncin.variables['lon'][:]
ncin.close()

3. Detrend

Here we used the detrend from scipy.signal. See help: signal.detrend?



In [4]:

    
nt,nlat,nlon = hgt.shape    
hgt = hgt.reshape((nt,nlat*nlon), order='F')
hgt_detrend = signal.detrend(hgt, axis=0, type='linear', bp=0)
hgt_detrend = hgt_detrend.reshape((nt,nlat,nlon), order='F')
print(hgt_detrend.shape)









    



(25L, 73L, 144L)

4. Carry out EOF analysis

4.1 Create an EOF solver to do the EOF analysis

Cosine of latitude weights are applied before the computation of EOFs.



In [5]:

    
wgts   = np.cos(np.deg2rad(lat))
wgts   = wgts.reshape(len(wgts), 1)
solver = Eof(hgt_detrend, weights=wgts)

4.3 Retrieve the leading EOFs

Expressed as the correlation between the leading PC time series and the input SST anomalies at each grid point, and the leading PC time series itself.



In [6]:

    
eof1 = solver.eofs(neofs=10)
pc1  = solver.pcs(npcs=10, pcscaling=0)
varfrac = solver.varianceFraction()
lambdas = solver.eigenvalues()

5. Visualize leading EOFs

5.1 Plot EOFs and PCs



In [7]:

    
parallels = np.arange(-90,90,30.)
meridians = np.arange(-180,180,30)

for i in range(0,2):
    fig = plt.figure(figsize=(12,9))
    plt.subplot(211)
    
    m = Basemap(projection='cyl', llcrnrlon=min(lon), llcrnrlat=min(lat), urcrnrlon=max(lon), urcrnrlat=max(lat))    
    x, y = m(*np.meshgrid(lon, lat))
    clevs = np.linspace(np.min(eof1[i,:,:].squeeze()), np.max(eof1[i,:,:].squeeze()), 21)
    cs = m.contourf(x, y, eof1[i,:,:].squeeze(), clevs, cmap=plt.cm.RdBu_r)
    m.drawcoastlines()  
    m.drawparallels(parallels, labels=[1,0,0,0])
    m.drawmeridians(meridians, labels=[1,0,0,1])

    cb = m.colorbar(cs, 'right', size='5%', pad='2%')
    cb.set_label('EOF', fontsize=12)
    plt.title('EOF ' + str(i+1), fontsize=16)

    plt.subplot(212)
    days = np.linspace(1979,2003,nt)
    plt.plot(days, pc1[:,i], linewidth=2)
    plt.axhline(0, color='k')
    plt.xlabel('Year')
    plt.ylabel('PC Amplitude')   
    plt.ylim(np.min(pc1.squeeze()), np.max(pc1.squeeze()))

5.2 Check variances explained by leading EOFs



In [8]:

    
plt.figure(figsize=(11,6))
eof_num = range(1, 16)
plt.plot(eof_num, varfrac[0:15], linewidth=2)
plt.plot(eof_num, varfrac[0:15], linestyle='None', marker="o", color='r', markersize=8)
plt.axhline(0, color='k')
plt.xticks(range(1, 16))
plt.title('Fraction of the total variance represented by each EOF')
plt.xlabel('EOF #')
plt.ylabel('Variance Fraction')
plt.xlim(1, 15)
plt.ylim(np.min(varfrac), np.max(varfrac)+0.01)









    Out[8]:





(1.448078e-16, 0.19716533482074738)

References

http://unidata.github.io/netcdf4-python/

John D. Hunter. Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 9, 90-95 (2007), DOI:10.1109/MCSE.2007.55

Stéfan van der Walt, S. Chris Colbert and Gaël Varoquaux. The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science & Engineering, 13, 22-30 (2011), DOI:10.1109/MCSE.2011.37

Kalnay et al.,The NCEP/NCAR 40-year reanalysis project, Bull. Amer. Meteor. Soc., 77, 437-470, 1996.



In [ ]: